Network Kernal Density Estimation (NKDE)

Published

April 10, 2024

Modified

April 10, 2024

1 Getting Started

pacman::p_load(arrow, dplyr, lubridate, ggplot2, tidyverse, sf, sfdep, sp, spNetwork, spdep, raster, spatstat, tmap, readxl, plotly)
osm_basemap <- tm_basemap(server = "OpenStreetMap.HOT")
imagery_basemap <- tm_basemap(server = "Esri.WorldImagery")

2. Importing Data

2.1 Asaptial Data

hk_census <- read_excel("data/aspatial/hkcensus.xlsx")

2.2 Geospatial Data

2.2.1 18 Districts in Hong Kong

district_18 <- st_read(dsn = "data/geospatial/hk_18Districts/",
                       layer = "HKDistrict18" )
Reading layer `HKDistrict18' from data source 
  `C:\Users\Kachel Lee\ka33rina\IS415-project\Analysis\data\geospatial\hk_18Districts' 
  using driver `ESRI Shapefile'
Simple feature collection with 18 features and 4 fields
Geometry type: MULTIPOLYGON
Dimension:     XY
Bounding box:  xmin: 12672060 ymin: 2529946 xmax: 12739620 ymax: 2579129
Projected CRS: WGS 84 / Pseudo-Mercator
sf_district_18 <- district_18 %>% st_transform(crs = 2326)
st_crs(sf_district_18)
Coordinate Reference System:
  User input: EPSG:2326 
  wkt:
PROJCRS["Hong Kong 1980 Grid System",
    BASEGEOGCRS["Hong Kong 1980",
        DATUM["Hong Kong 1980",
            ELLIPSOID["International 1924",6378388,297,
                LENGTHUNIT["metre",1]]],
        PRIMEM["Greenwich",0,
            ANGLEUNIT["degree",0.0174532925199433]],
        ID["EPSG",4611]],
    CONVERSION["Hong Kong 1980 Grid",
        METHOD["Transverse Mercator",
            ID["EPSG",9807]],
        PARAMETER["Latitude of natural origin",22.3121333333333,
            ANGLEUNIT["degree",0.0174532925199433],
            ID["EPSG",8801]],
        PARAMETER["Longitude of natural origin",114.178555555556,
            ANGLEUNIT["degree",0.0174532925199433],
            ID["EPSG",8802]],
        PARAMETER["Scale factor at natural origin",1,
            SCALEUNIT["unity",1],
            ID["EPSG",8805]],
        PARAMETER["False easting",836694.05,
            LENGTHUNIT["metre",1],
            ID["EPSG",8806]],
        PARAMETER["False northing",819069.8,
            LENGTHUNIT["metre",1],
            ID["EPSG",8807]]],
    CS[Cartesian,2],
        AXIS["northing (N)",north,
            ORDER[1],
            LENGTHUNIT["metre",1]],
        AXIS["easting (E)",east,
            ORDER[2],
            LENGTHUNIT["metre",1]],
    USAGE[
        SCOPE["Cadastre, engineering survey, topographic mapping (large scale)."],
        AREA["China - Hong Kong - onshore and offshore."],
        BBOX[22.13,113.76,22.58,114.51]],
    ID["EPSG",2326]]
tmap_mode("plot")
d18_map <- tm_shape(sf_district_18) +
  tm_fill(col = "ENAME", title = "District", legend.show = FALSE) +
  tm_borders(col = "black", lwd = 0.5) +
  tm_text("ID", size = 0.5, col = "black") +
  tm_layout(frame = FALSE)

d18_map

2.2.2 Recycling Points in Hong Kong

cp <- read_csv("data/aspatial/hkrecyclepoints.csv")

Change the geographic location

cp_sf <- st_as_sf(cp, 
                  coords = c("lgt","lat"), 
                  crs = 4326) %>% 
  st_transform(crs= 2326)
summary(cp_sf)
     cp_id       district_id         address_en        address2_en       
 Min.   :  283   Length:6520        Length:6520        Length:6520       
 1st Qu.: 3615   Class :character   Class :character   Class :character  
 Median : 6652   Mode  :character   Mode  :character   Mode  :character  
 Mean   : 6127                                                           
 3rd Qu.: 8516                                                           
 Max.   :10764                                                           
  waste_type           legend                   geometry   
 Length:6520        Length:6520        POINT        :6520  
 Class :character   Class :character   epsg:2326    :   0  
 Mode  :character   Mode  :character   +proj=tmer...:   0  
                                                           
                                                           
                                                           
cp_sf_1 <- cp_sf %>%
  mutate(district_id = toupper(str_replace_all(district_id, "_", " ")))
recycling_bins <- subset(cp_sf_1, legend == "Recycling Bins at Public Place")
recycling_spots <- subset(cp_sf_1, legend == "Recycling Spots")
private_collection_points <- subset(cp_sf_1, legend == "Private Collection Points (e.g. housing estates, shopping centres)")
ngo_collection_points <- subset(cp_sf_1, legend == "NGO Collection Points")
recycling_stations <- subset(cp_sf_1, legend == "Recycling Stations/Recycling Stores")
street_corner_recycling_shops <- subset(cp_sf_1, legend == "Street Corner Recycling Shops")
smart_bins <- subset(cp_sf_1, legend == "Smart Bin")
recycling_spots_cp <- st_join(sf_district_18, recycling_spots)
ngo_cp <- st_join(sf_district_18, ngo_collection_points)
pcp_joined_data <- st_join(sf_district_18, private_collection_points)
recycling_bins_cp <- st_join(sf_district_18, recycling_bins)
recycling_stations_cp <- st_join(sf_district_18, recycling_stations)
street_corner_cp <- st_join(sf_district_18, street_corner_recycling_shops)
smart_bins_cp <- st_join(sf_district_18, smart_bins)
private_collection_points_by_district <- pcp_joined_data %>%
  group_by(ENAME) %>%
  summarize(total_pcp = n())

ngo_cp_by_district <- ngo_cp %>%
  group_by(ENAME) %>%
  summarize(total_ngo_cp = n())

recycling_spots_by_district <- recycling_spots_cp %>%
  group_by(ENAME) %>%
  summarize(total_recycling_spots = n())

recycling_bins_by_district <- recycling_bins_cp %>%
  group_by(ENAME) %>%
  summarize(total_recycling_bins = n())

recycling_stations_by_district <- recycling_stations_cp %>%
  group_by(ENAME) %>%
  summarize(total_recycling_stations = n())

street_corner_shops_by_district <- street_corner_cp %>%
  group_by(ENAME) %>%
  summarize(total_street_corner = n())

smart_bins_by_district <- smart_bins_cp %>%
  group_by(ENAME) %>%
  summarize(total_smart_bins = n())

2.3 Road Data in Hong Kong

road_data <- st_read(dsn = "data/geospatial/china-latest-free.shp", 
                     layer = "gis_osm_roads_free_1")

Transform the Data into Hong Kong Projection System

road_data_2326 <- st_transform(road_data, 2326)
roads_in_hk <- st_intersection(road_data_2326, sf_district_18)
write_rds(roads_in_hk, "data/rds/sf_roads_in_hk.rds")
roads_in_hk <- read_rds("data/rds/sf_roads_in_hk.rds")

3. Visualizing Recycling Points in Hong Kong

pcp_map <- tm_shape(private_collection_points_by_district) +
  tm_fill(col = "total_pcp") +
  tm_borders() +
  tm_layout(legend.show = TRUE, main.title = "Distribution of Private Collection Points by District",main.title.position = "center", main.title.size = 0.75)+
  tm_scale_bar()
ngo_cp_map <- tm_shape(ngo_cp_by_district) +
  tm_fill(col = "total_ngo_cp") +
  tm_borders() +
  tm_layout(legend.show = TRUE, main.title = "Distribution of NGO Points by District", main.title.position = "center", main.title.size = 0.75)+
  tm_scale_bar()
recycling_spots_map <- tm_shape(recycling_spots_by_district) +
  tm_fill(col = "total_recycling_spots") +
  tm_borders() +
  tm_layout(legend.show = TRUE, main.title = "Distribution of Recycling Spots by District",main.title.position = "center", main.title.size = 0.75)+
  tm_scale_bar()
recycling_bins_map <- tm_shape(recycling_bins_by_district) +
  tm_fill(col = "total_recycling_bins") +
  tm_borders() +
  tm_layout(legend.show = TRUE, main.title = "Distribution of Recycling Bins by District", main.title.position = "center", main.title.size = 0.75)+
  tm_scale_bar()
recycling_stations_map <- tm_shape(recycling_stations_by_district) +
  tm_fill(col = "total_recycling_stations") +
  tm_borders() +
  tm_layout(legend.show = TRUE, main.title = "Distribution of Recycling Stations by District", main.title.position = "center", main.title.size = 0.75)+
  tm_scale_bar()
street_corner_shops_map <- tm_shape(street_corner_shops_by_district) +
  tm_fill(col = "total_street_corner") +
  tm_borders() +
  tm_layout(legend.show = TRUE,main.title = "Distribution of Street Corner Shops by District", main.title.position = "center", main.title.size = 0.75)+
  tm_scale_bar()
smart_bin_map <- tm_shape(smart_bins_by_district) +
  tm_fill(col = "total_smart_bins") +
  tm_borders() +
  tm_layout(legend.show = TRUE,main.title = "Distribution of Smart Bins by District",main.title.position = "center", main.title.size = 0.75)+
  tm_scale_bar()
tmap_arrange(pcp_map,recycling_spots_map, recycling_bins_map,ngo_cp_map, recycling_stations_map, street_corner_shops_map, smart_bin_map, ncol=2, nrow = 4)

Since we want to conduct a more accurate analysis on the recycling points, We will focus on the recycling facilities type with a higher number of data - Private Collection Points and Recycling Bins.

3.1 Private Collection Points

wm_q_pcp <- private_collection_points_by_district %>%
  mutate(nb = st_contiguity(geometry, queen = TRUE),
         wt = st_weights(nb,
                         style = "W",
                         allow_zero = TRUE),
         .before = 1) 
wm_q_pcp
Simple feature collection with 18 features and 4 fields
Geometry type: MULTIPOLYGON
Dimension:     XY
Bounding box:  xmin: 801026 ymin: 801678 xmax: 863539.9 ymax: 846903.1
Projected CRS: Hong Kong 1980 Grid System
# A tibble: 18 × 5
   nb        wt        ENAME             total_pcp                      geometry
 * <nb>      <list>    <chr>                 <int>            <MULTIPOLYGON [m]>
 1 <int [2]> <dbl [2]> CENTRAL & WESTERN       414 (((833047.3 816838.9, 833488…
 2 <int [2]> <dbl [2]> EASTERN                 200 (((843535.7 812736.3, 843530…
 3 <int [1]> <dbl [1]> ISLANDS                  46 (((810029 801683.4, 810019.4…
 4 <int [5]> <dbl [5]> KOWLOON CITY            243 (((836281 823326.3, 836283 8…
 5 <int [3]> <dbl [3]> KWAI TSING               87 (((827834 820693.1, 827834 8…
 6 <int [3]> <dbl [3]> KWUN TONG               119 (((843155.9 816368.9, 843154…
 7 <int [2]> <dbl [2]> NORTH                    97 (((852615.5 841161.9, 852615…
 8 <int [4]> <dbl [4]> SAI KUNG                163 (((840825.4 823785.4, 840827…
 9 <int [7]> <dbl [7]> SHA TIN                 194 (((839924.4 832031.6, 839924…
10 <int [4]> <dbl [4]> SHAM SHUI PO            190 (((834919.6 823272.4, 834921…
11 <int [3]> <dbl [3]> SOUTHERN                181 (((830176.1 815000.2, 830255…
12 <int [5]> <dbl [5]> TAI PO                  121 (((846104.6 830686.6, 846265…
13 <int [6]> <dbl [6]> TSUEN WAN                92 (((821762.8 819430, 821765.7…
14 <int [2]> <dbl [2]> TUEN MUN                165 (((811619.1 831909.9, 811639…
15 <int [3]> <dbl [3]> WAN CHAI                198 (((838663.3 815002.8, 838700…
16 <int [4]> <dbl [4]> WONG TAI SIN             82 (((836530.5 823327.5, 836533…
17 <int [2]> <dbl [2]> YAU TSIM MONG           116 (((835145.2 817859.4, 835137…
18 <int [4]> <dbl [4]> YUEN LONG               178 (((811708.8 831974.4, 811719…
set.seed(1234)
global_moran_perm(wm_q_pcp$total_pcp,
                       wm_q_pcp$nb,
                       wm_q_pcp$wt,
                  zero.policy = TRUE,
                  nsim = 999)

    Monte-Carlo simulation of Moran I

data:  x 
weights: listw  
number of simulations + 1: 1000 

statistic = 0.16065, observed rank = 934, p-value = 0.132
alternative hypothesis: two.sided
lisa_pcp <- wm_q_pcp %>% 
  mutate(local_moran = local_moran(
    total_pcp, nb, wt, zero.policy = TRUE, nsim = 99),
         .before = 1) %>%
  unnest(local_moran)
lisa_sig_pcp <- lisa_pcp  %>%
  filter(p_ii < 0.05)

tmap_mode("plot")
tm_shape(lisa_pcp) +
  tm_polygons() +
  tm_borders(alpha = 0.7) +
tm_shape(lisa_pcp) +
  tm_fill("mean") + 
  tm_borders(alpha = 0.4)

3.2 Recycling Bins

wm_q_rbins <- recycling_bins_by_district %>%
  mutate(nb = st_contiguity(geometry, queen = TRUE),
         wt = st_weights(nb,
                         style = "W",
                         allow_zero = TRUE),
         .before = 1) 
wm_q_rbins
Simple feature collection with 18 features and 4 fields
Geometry type: MULTIPOLYGON
Dimension:     XY
Bounding box:  xmin: 801026 ymin: 801678 xmax: 863539.9 ymax: 846903.1
Projected CRS: Hong Kong 1980 Grid System
# A tibble: 18 × 5
   nb        wt        ENAME      total_recycling_bins                  geometry
 * <nb>      <list>    <chr>                     <int>        <MULTIPOLYGON [m]>
 1 <int [2]> <dbl [2]> CENTRAL &…                   93 (((833047.3 816838.9, 83…
 2 <int [2]> <dbl [2]> EASTERN                     177 (((843535.7 812736.3, 84…
 3 <int [1]> <dbl [1]> ISLANDS                     219 (((810029 801683.4, 8100…
 4 <int [5]> <dbl [5]> KOWLOON C…                   92 (((836281 823326.3, 8362…
 5 <int [3]> <dbl [3]> KWAI TSING                  130 (((827834 820693.1, 8278…
 6 <int [3]> <dbl [3]> KWUN TONG                   133 (((843155.9 816368.9, 84…
 7 <int [2]> <dbl [2]> NORTH                       286 (((852615.5 841161.9, 85…
 8 <int [4]> <dbl [4]> SAI KUNG                    302 (((840825.4 823785.4, 84…
 9 <int [7]> <dbl [7]> SHA TIN                     214 (((839924.4 832031.6, 83…
10 <int [4]> <dbl [4]> SHAM SHUI…                  110 (((834919.6 823272.4, 83…
11 <int [3]> <dbl [3]> SOUTHERN                    115 (((830176.1 815000.2, 83…
12 <int [5]> <dbl [5]> TAI PO                      282 (((846104.6 830686.6, 84…
13 <int [6]> <dbl [6]> TSUEN WAN                   193 (((821762.8 819430, 8217…
14 <int [2]> <dbl [2]> TUEN MUN                    195 (((811619.1 831909.9, 81…
15 <int [3]> <dbl [3]> WAN CHAI                    121 (((838663.3 815002.8, 83…
16 <int [4]> <dbl [4]> WONG TAI …                   60 (((836530.5 823327.5, 83…
17 <int [2]> <dbl [2]> YAU TSIM …                  147 (((835145.2 817859.4, 83…
18 <int [4]> <dbl [4]> YUEN LONG                   327 (((811708.8 831974.4, 81…
set.seed(1234)
global_moran_perm(wm_q_rbins$total_recycling_bins,
                       wm_q_rbins$nb,
                       wm_q_rbins$wt,
                  zero.policy = TRUE,
                  nsim = 999)

    Monte-Carlo simulation of Moran I

data:  x 
weights: listw  
number of simulations + 1: 1000 

statistic = 0.48951, observed rank = 998, p-value = 0.004
alternative hypothesis: two.sided
lisa_rbins <- wm_q_rbins %>% 
  mutate(local_moran = local_moran(
    total_recycling_bins, nb, wt, zero.policy = TRUE, nsim = 99),
         .before = 1) %>%
  unnest(local_moran)
lisa_sig_rbins <- lisa_rbins  %>%
  filter(p_ii < 0.05)

tmap_mode("plot")
tm_shape(lisa_rbins) +
  tm_polygons() +
  tm_borders(alpha = 0.7) +
tm_shape(lisa_rbins) +
  tm_fill("mean") + 
  tm_borders(alpha = 0.4)

Considering the Hot-Spot Analysis on the Private Collection Points and Public Recycling bins, they are maininly locate far from the residential area. Therefore, we decide to refer to the census data in Hong Kong.

ggplot(hk_census, aes(x = DNAME, y = POPULATION)) +
  geom_bar(stat = "identity", fill = "steelblue", color = "black") +
  geom_text(aes(label = POPULATION), vjust = -0.5) +
  theme_minimal() +
  labs(title = "Population by District",
       x = "District Name",
       y = "Population") +
  theme(axis.text.x = element_text(angle = 45, hjust = 1))

Since Sha Tin has the most population, we will focus on Sha Tin district. Moreover, we want to choose another district with both commercial and residential area. Therefore, we will be analyzing on Yau Tsim Mong Area as well.

sha_tin_recycling_points <- cp_sf_1 %>%
  filter(district_id == "SHA TIN" & 
         (legend == "Private Collection Points (e.g. housing estates, shopping centres)" | legend == "Recycling Bins at Public Place"))
ytm_recycling_points <- cp_sf_1 %>%
  filter(district_id == "YAU TSIM MONG" & 
         (legend == "Private Collection Points (e.g. housing estates, shopping centres)" | legend == "Recycling Bins at Public Place"))
wanchai_recycling_points <- cp_sf_1 %>%
  filter(district_id == "WAN CHAI" & 
         (legend == "Private Collection Points (e.g. housing estates, shopping centres)" | legend == "Recycling Bins at Public Place"))

4. Sha Tin District

shatin <- sf_district_18 %>%
  filter(ENAME == "SHA TIN")
plot(st_geometry((shatin)))

roads_in_shatin <- st_intersection(roads_in_hk, shatin)
write_rds(roads_in_shatin, "data/rds/roads_in_shatin.rds")
roads_in_shatin <- read_rds("data/rds/roads_in_shatin.rds")
private_cp_shatin <- st_intersection(private_collection_points, shatin)
public_cp_shatin <- st_intersection(recycling_bins, shatin)
cp_shatin <- st_intersection(sha_tin_recycling_points, shatin)
write_rds(private_cp_shatin, "data/rds/private_cp_shatin.rds")
write_rds(public_cp_shatin, "data/rds/public_cp_shatin.rds")
write_rds(cp_shatin, "data/rds/cp_shatin.rds")
private_cp_shatin <- read_rds("data/rds/private_cp_shatin.rds")
public_cp_shatin <- read_rds("data/rds/public_cp_shatin.rds")
cp_shatin <- read_rds("data/rds/cp_shatin.rds")
plot(st_geometry(private_cp_shatin))

plot(st_geometry(public_cp_shatin))

plot(st_geometry(cp_shatin))

unique_types <- unique(st_geometry_type(roads_in_shatin))
unique_types
[1] POINT              GEOMETRYCOLLECTION LINESTRING         MULTILINESTRING   
[5] MULTIPOINT        
18 Levels: GEOMETRY POINT LINESTRING POLYGON MULTIPOINT ... TRIANGLE
roads_in_shatin <- roads_in_shatin %>%
  filter(st_geometry_type(geometry) %in% c("LINESTRING", "MULTILINESTRING"))
if ("LINESTRING" %in% unique_types) {
  roads_in_shatin <- st_cast(roads_in_shatin, "LINESTRING")
} else {
  # handle the case when no linestrings are found
  stop("No linestrings found in roads_in_shatin")
}
tmap_mode('plot')
tm_shape(roads_in_shatin, geometry_type = "lines") + 
  tm_lines()

tmap_mode('plot')

# Plot roads
tm_shape(roads_in_shatin, geometry_type = "lines") + 
  tm_lines(lwd = 1, col = "blue") +
  
  # Plot Punggol area boundary
  tm_shape(shatin) +
  tm_borders() +
  
  # Plot Grab origin data
  tm_shape(private_cp_shatin) +
  tm_dots()

tmap_mode("view")

# Create the map using tmap
cp_shatin_map <- osm_basemap + 
  tm_shape(cp_shatin) +
  tm_dots(col = "legend", palette = c("blue", "red"), border.col = "black", size = 0.05) +
  tm_layout(title = "Recycling Bin Locations in Sha Tin")

cp_shatin_map
roads_lines_shatin <- roads_in_shatin[st_geometry_type(roads_in_shatin) == "LINESTRING", ]

# Apply lixelize_lines with mindist
lixels_shatin <- lixelize_lines(roads_lines_shatin,5000, mindist = 2500)
samples_shatin <- lines_center(roads_lines_shatin)

4.1 Private Collection Points in Sha Tin

densities_shatin_private <- nkde(roads_lines_shatin, 
                  events = private_cp_shatin,
                  w = rep(1,nrow(private_cp_shatin)),
                  samples = samples_shatin,
                  kernel_name = "quartic",
                  bw = 300, 
                  div= "bw", 
                  method = "simple", 
                  digits = 1, 
                  tol = 1,
                  grid_shape = c(1,1), 
                  max_depth = 8,
                  agg = 5, #we aggregate events within a 5m radius (faster calculation)
                  sparse = TRUE,
                  verbose = FALSE)
samples_shatin$density_private <- densities_shatin_private
lixels_shatin$density_private <- densities_shatin_private
# rescaling to help the mapping
samples_shatin$density_private <- samples_shatin$density_private*1000
lixels_shatin$density_private <- lixels_shatin$density_private*1000
tmap_mode('view')

shatin_density_private <- osm_basemap+
tm_shape(lixels_shatin)+
  tm_lines(col="density_private")+
tm_shape(private_cp_shatin)+
  tm_dots()

shatin_density_private

4.2 Public Recycling Bins in Sha Tin

densities_shatin_public <- nkde(roads_lines_shatin, 
                  events = public_cp_shatin,
                  w = rep(1,nrow(public_cp_shatin)),
                  samples = samples_shatin,
                  kernel_name = "quartic",
                  bw = 300, 
                  div= "bw", 
                  method = "simple", 
                  digits = 1, 
                  tol = 1,
                  grid_shape = c(1,1), 
                  max_depth = 8,
                  agg = 5, #we aggregate events within a 5m radius (faster calculation)
                  sparse = TRUE,
                  verbose = FALSE)
samples_shatin$density_public <- densities_shatin_public
lixels_shatin$density_public <- densities_shatin_public
# rescaling to help the mapping
samples_shatin$density_public <- samples_shatin$density_public*1000
lixels_shatin$density_public <- lixels_shatin$density_public*1000
tmap_mode('view')


shatin_density_public <- osm_basemap+
tm_shape(lixels_shatin)+
  tm_lines(col="density_public")+
tm_shape(public_cp_shatin)+
  tm_dots()

shatin_density_public

4.3 Both Private and Public Bins in Sha Tin Districts

densities_shatin_ppcp <- nkde(roads_lines_shatin, 
                  events = cp_shatin,
                  w = rep(1,nrow(cp_shatin)),
                  samples = samples_shatin,
                  kernel_name = "quartic",
                  bw = 300, 
                  div= "bw", 
                  method = "simple", 
                  digits = 1, 
                  tol = 1,
                  grid_shape = c(1,1), 
                  max_depth = 8,
                  agg = 5, #we aggregate events within a 5m radius (faster calculation)
                  sparse = TRUE,
                  verbose = FALSE)
samples_shatin$density_ppcp <- densities_shatin_ppcp
lixels_shatin$density_ppcp <- densities_shatin_ppcp
# rescaling to help the mapping
samples_shatin$density_ppcp <- samples_shatin$density_ppcp*1000
lixels_shatin$density_ppcp <- lixels_shatin$density_ppcp*1000
tmap_mode('view')

shatin_density_ppcp <- osm_basemap +
  tm_shape(lixels_shatin) +
  tm_lines(col = "density_ppcp") +
  tm_shape(cp_shatin) +
  tm_layout(
    legend.position = c("left", "top"),
    legend.text.size = 0.5,  # Adjust this value as needed for smaller legend text
    legend.title.size = 0.6  # Adjust this value as needed for smaller legend title
  )+
  tm_dots(col = "legend", palette = c("blue", "red"), border.col = "black", size = 0.05) 

shatin_density_ppcp

5. Yau Tsim Mong Area

ytm <- sf_district_18 %>%
  filter(ENAME == "YAU TSIM MONG")
plot(st_geometry((ytm)))

roads_in_ytm <- st_intersection(roads_in_hk, ytm)
write_rds(roads_in_ytm, "data/rds/roads_in_ytm.rds")
roads_in_ytm <- read_rds("data/rds/roads_in_ytm.rds")
private_cp_ytm <- st_intersection(private_collection_points, ytm)
public_cp_ytm <- st_intersection(recycling_bins, ytm)
cp_ytm <- st_intersection(ytm_recycling_points, ytm)
write_rds(private_cp_ytm, "data/rds/private_cp_ytm.rds")
write_rds(public_cp_ytm, "data/rds/public_cp_ytm.rds")
write_rds(cp_ytm, "data/rds/cp_ytm.rds")
private_cp_ytm <- read_rds("data/rds/private_cp_ytm.rds")
public_cp_ytm <- read_rds("data/rds/public_cp_ytm.rds")
cp_ytm <-  read_rds("data/rds/cp_ytm.rds")
plot(st_geometry(private_cp_ytm))

plot(st_geometry(public_cp_ytm))

unique_types <- unique(st_geometry_type(roads_in_ytm))
unique_types
[1] LINESTRING         MULTILINESTRING    POINT              GEOMETRYCOLLECTION
18 Levels: GEOMETRY POINT LINESTRING POLYGON MULTIPOINT ... TRIANGLE
roads_in_ytm <- roads_in_ytm %>%
  filter(st_geometry_type(geometry) %in% c("LINESTRING", "MULTILINESTRING"))
if ("LINESTRING" %in% unique_types) {
  roads_in_ytm <- st_cast(roads_in_ytm, "LINESTRING")
} else {
  # handle the case when no linestrings are found
  stop("No linestrings found in roads_in_ytm")
}
tmap_mode('plot')
tm_shape(roads_in_ytm, geometry_type = "lines") + 
  tm_lines()

5.1 Private Collection Poitns in YTM

tmap_mode('plot')

# Plot roads
tm_shape(roads_in_ytm, geometry_type = "lines") + 
  tm_lines(lwd = 1, col = "blue") +
  
  tm_shape(ytm) +
  tm_borders() +

  tm_shape(private_cp_ytm) +
  tm_dots()
tmap_mode("view")

# Create the map using tmap
cp_ytm_map <- osm_basemap + 
  tm_shape(cp_ytm) +
  tm_dots(col = "legend", palette = c("blue", "red"), border.col = "black", size = 0.05) +
  tm_layout(title = "Recycling Bin Locations in Yau Tsim Mong")

cp_shatin_map
roads_lines_ytm <- roads_in_ytm[st_geometry_type(roads_in_ytm) == "LINESTRING", ]

# Apply lixelize_lines with mindist
lixels_ytm <- lixelize_lines(roads_lines_ytm,5000, mindist = 2500)
samples_ytm <- lines_center(lixels_ytm)
densities_ytm_private <- nkde(roads_lines_ytm, 
                  events = private_cp_ytm,
                  w = rep(1,nrow(private_cp_ytm)),
                  samples = samples_ytm,
                  kernel_name = "quartic",
                  bw = 300, 
                  div= "bw", 
                  method = "simple", 
                  digits = 1, 
                  tol = 1,
                  grid_shape = c(1,1), 
                  max_depth = 8,
                  agg = 5, #we aggregate events within a 5m radius (faster calculation)
                  sparse = TRUE,
                  verbose = FALSE)
samples_ytm$density_private <- densities_ytm_private
lixels_ytm$density_private <- densities_ytm_private
# rescaling to help the mapping
samples_ytm$density_private <- samples_ytm$density_private*1000
lixels_ytm$density_private <- lixels_ytm$density_private*1000
tmap_mode('view')

ytm_density_private <- osm_basemap + 
  tm_shape(lixels_ytm)+
  tm_lines(col="density_private")+
tm_shape(private_cp_ytm)+
  tm_dots()

ytm_density_private

5.2 Public Recycling Bins in YTM

densities_ytm_public <- nkde(roads_lines_ytm, 
                  events = public_cp_ytm,
                  w = rep(1,nrow(public_cp_ytm)),
                  samples = samples_ytm,
                  kernel_name = "quartic",
                  bw = 300, 
                  div= "bw", 
                  method = "simple", 
                  digits = 1, 
                  tol = 1,
                  grid_shape = c(1,1), 
                  max_depth = 8,
                  agg = 5, #we aggregate events within a 5m radius (faster calculation)
                  sparse = TRUE,
                  verbose = FALSE)
samples_ytm$density_public <- densities_ytm_public
lixels_ytm$density_public <- densities_ytm_public
# rescaling to help the mapping
samples_ytm$density_public <- samples_ytm$density_public*1000
lixels_ytm$density_public <- lixels_ytm$density_public*1000
tmap_mode('view')

ytm_density_public <-osm_basemap + 
tm_shape(lixels_ytm)+
  tm_lines(col="density_public")+
tm_shape(private_cp_ytm)+
  tm_dots()

ytm_density_public

5.3 Both Public and Private Recycling Bins in YTM

densities_ytm_ppcp <- nkde(roads_lines_ytm, 
                  events = cp_ytm,
                  w = rep(1,nrow(cp_ytm)),
                  samples = samples_ytm,
                  kernel_name = "quartic",
                  bw = 300, 
                  div= "bw", 
                  method = "simple", 
                  digits = 1, 
                  tol = 1,
                  grid_shape = c(1,1), 
                  max_depth = 8,
                  agg = 5, #we aggregate events within a 5m radius (faster calculation)
                  sparse = TRUE,
                  verbose = FALSE)
samples_ytm$density_ppcp <- densities_ytm_ppcp
lixels_ytm$density_ppcp <- densities_ytm_ppcp
# rescaling to help the mapping
samples_ytm$density_ppcp <- samples_ytm$density_ppcp*1000
lixels_ytm$density_ppcp <- lixels_ytm$density_ppcp*1000
tmap_mode('view')

ytm_density_ppcp <-osm_basemap + 
  tm_shape(lixels_ytm)+
  tm_lines(col="density_ppcp")+
  tm_shape(cp_ytm)+
  tm_dots(col = "legend", palette = c("blue", "red"), border.col = "black", size = 0.05)
  

ytm_density_ppcp

6. Wan Chai (Hong Kong Island)

wanchai <- sf_district_18 %>%
  filter(ENAME == "WAN CHAI")
plot(st_geometry((wanchai)))

roads_in_wanchai <- st_intersection(roads_in_hk, wanchai)
write_rds(roads_in_wanchai, "data/rds/roads_in_wanchai.rds")
roads_in_wanchai <- read_rds("data/rds/roads_in_wanchai.rds")
private_cp_wanchai <- st_intersection(private_collection_points, wanchai)
public_cp_wanchai<- st_intersection(recycling_bins, wanchai)
cp_wanchai <- st_intersection(wanchai_recycling_points, wanchai)
write_rds(private_cp_wanchai, "data/rds/private_cp_wanchai.rds")
write_rds(public_cp_wanchai, "data/rds/public_cp_wanchai.rds")
write_rds(cp_wanchai, "data/rds/cp_wanchai.rds")
private_cp_wanchai <- read_rds("data/rds/private_cp_wanchai.rds")
public_cp_wanchai <- read_rds("data/rds/public_cp_wanchai.rds")
cp_wanchai <- read_rds("data/rds/cp_wanchai.rds")
plot(st_geometry(private_cp_wanchai))

plot(st_geometry(cp_wanchai))

unique_types <- unique(st_geometry_type(roads_in_wanchai))
unique_types
[1] GEOMETRYCOLLECTION POINT              LINESTRING         MULTIPOINT        
[5] MULTILINESTRING   
18 Levels: GEOMETRY POINT LINESTRING POLYGON MULTIPOINT ... TRIANGLE
roads_in_wanchai <- roads_in_wanchai %>%
  filter(st_geometry_type(geometry) %in% c("LINESTRING", "MULTILINESTRING"))
if ("LINESTRING" %in% unique_types) {
  roads_in_wanchai <- st_cast(roads_in_wanchai, "LINESTRING")
} else {
  # handle the case when no linestrings are found
  stop("No linestrings found in roads_in_ytm")
}
tmap_mode('plot')
tm_shape(roads_in_wanchai, geometry_type = "lines") + 
  tm_lines()

6.3 Both Private and Public Bins in Wan Chai Districts

roads_lines_wanchai <- roads_in_wanchai[st_geometry_type(roads_in_wanchai) == "LINESTRING", ]

# Apply lixelize_lines with mindist
lixels_wanchai <- lixelize_lines(roads_lines_wanchai,5000, mindist = 2500)
samples_wanchai <- lines_center(lixels_wanchai)
densities_wanchai_ppcp <- nkde(roads_lines_wanchai, 
                  events = cp_wanchai,
                  w = rep(1,nrow(cp_wanchai)),
                  samples = samples_wanchai,
                  kernel_name = "quartic",
                  bw = 300, 
                  div= "bw", 
                  method = "simple", 
                  digits = 1, 
                  tol = 1,
                  grid_shape = c(1,1), 
                  max_depth = 8,
                  agg = 5, #we aggregate events within a 5m radius (faster calculation)
                  sparse = TRUE,
                  verbose = FALSE)
samples_wanchai$density_ppcp <- densities_wanchai_ppcp
lixels_wanchai$density_ppcp <- densities_wanchai_ppcp
# rescaling to help the mapping
samples_wanchai$density_ppcp <- samples_wanchai$density_ppcp*1000
lixels_wanchai$density_ppcp <- lixels_wanchai$density_ppcp*1000
tmap_mode('view')

wanchai_density_ppcp <- osm_basemap +
  tm_shape(lixels_wanchai) +
  tm_lines(col = "density_ppcp") +
  tm_shape(cp_wanchai) +
  tm_layout(
    legend.position = c("left", "top"),
    legend.text.size = 0.5,  # Adjust this value as needed for smaller legend text
    legend.title.size = 0.6  # Adjust this value as needed for smaller legend title
  )+
  tm_dots(col = "legend", palette = c("blue", "red"), border.col = "black", size = 0.05) 

wanchai_density_ppcp